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The Polarizable Continuum Model (PCM) can be used in conjunction with Density Functional Theory (DFT) 
and its time-dependent extension (TDDFT) to simulate the electronic and optical properties of molecules 
and nanoparticles immersed in a dielectric environment, typically liquid solvents. In this contribution, we 
develop a methodology to account for solvation effects in real-space (and real-time) (TD)DFT calculations. 
The boundary elements method is used to calculate the solvent reaction potential in terms of the apparent 
charges that spread over the Van der Waals solute surface. In a real-space representation this potential 
may exhibit a Coulomb singularity at grid points that are close to the cavity surface. We propose a simple 
approach to regularize such singularity by using a set of spherical Gaussian functions to distribute the apparent 
charges. We have implemented the proposed method in the Octopus code and present results for the 
electrostatic contribution to the solvation free energies and solvatochromic shifts for a representative set of 
organic molecules in water. 


I. INTRODUCTION 

Most of the innovative applications investigated in 
nanosciences rely on the physicochemical properties of 
novel materials which can be modulated in a non-trivial 
manner due to the interaction with surrounding environ¬ 
ments like solvent solutions, mesoscopic nanostructures 
and surfaces. The necessity of considering dielectric ef¬ 
fects on the opto-electronic and dynamical properties of 
molecular systems is exemplified by a vast scientific liter¬ 
ature addressing applications in different fields such as 
dye-sensitized solar cells^ water-soluble functionalized 
fullerens for biotechnology^^ porphyrin-based complexes 
used as artificial photosynthetic reaction center- 4 or as 
efficient sensitizers for photodynamic therapies^ among 
others. 

A widely used approach to describe molecules in con¬ 
densed phase consists in the so called focussed models 
where the system is divided into parts which are treated 
at different levels of accuracy^ An important subset 
within these models is constituted by the different for¬ 
mulations of the Polarizable Continuum Model (PCM) 7 
comprising the original dielectric PCM (D-PCM)&2, the 
conductor-like screening model (C-PCM) 10,1 - and the in¬ 
tegral equation formalism (IEF-PCM) j 12 i 13 All of them 
rely on three basic points as illustrated in Fig. [TJ i) 
the solute molecule is treated at the quantum mechanical 
level, ii) it is hosted by a cavity that exclude the solvent 
and contains the largest possible amount of the solute 
charge density and iii) the solvent is considered a contin¬ 
uous polarizable medium characterized by its frequency- 
dependent dielectric function e(w). 

PCM implementations are available among the most 
popular Quantum Chemistry programs using basis sets 
such as plane waves and atomic orbitals.— ‘ 14r — However, 
in the last years, with the rapid increase of computational 
power, real-space finite-difference methods have gained 


a lot of attention to perform first-principles electronic 
structure calculations*- 1 ™ - — They are numerically robust 
and very accurate since molecular orbitals defined on a 
spatial grid have the largest flexibility to take the proper 
values in both the intra- and the inter-atomic regions. 
Moreover, the potentials can be directly calculated in 
real space without the need of using basis sets, and arbi¬ 
trary boundary conditions can be used to simulate actual 
experimental environments. 

In this paper, we present a methodology to account for 
solvation effects within real-space and real-time calcula¬ 
tions and demonstrate that solvation energies and solva¬ 
tochromic shifts can be calculated to chemical accuracy 
as compared with state-of-the-art quantum-chemical ba¬ 
sis sets methods. Furthermore, we show that this method 
allows to investigate solvent effects in the real-time elec¬ 
tron dynamics of photoexcited states. In particular, we 
propose a simple method to regularize the Coulomb sin¬ 
gularity in the solvent reaction potential that may arise 
for grid points in the simulation domain that are infinites¬ 
imally close to the discretized solute-solvent interface. To 
that aim, we introduced a set of normalized spherical 
Gaussian functions to smooth the discretized polariza¬ 
tion charge density. This idea is very similar in spirit 
to the continuous surface charge formalism described by 
Scalmani et al. in Ref. [l5| even though the motivations 
are different. Within their approach the discretized PCM 
integral equations have been re-written to avoid discon¬ 
tinuities in the calculated potential energy surfaces for 
molecules in solution while in the present work we have 
focused on deriving a simple regularized expression for 
the solvent reaction potential. The method is imple¬ 
mented in the GPL license Octopus code— 

We use the regularized PCM potential to calculate the 
electrostatic contribution to the solvation free energies 
and the optical absorption response for a set of sim¬ 
ple organic molecules in water. To calculate the ground 
state energy of the solute systems, we use Density Func- 
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continuum solvent e(o)) 


FIG. 1. Schematic representation of the Nitrobenzene 
molecule embedded in a continuum solvent with a frequency- 
dependent dielectric constant e(w). V denotes the surface 
of the cavity hosting the solute system which is treated 
quantum-mechanically. 


tional Theory (DFT) by solving in real-space the self- 
consistent Kohn-Sham (KS) equations. The optical ab¬ 
sorption spectra is calculated by exciting the system with 
a sudden external perturbation and propagate the KS 
states in real-timei 21 ’ 22 In this work the real-time sol¬ 
vent response is calculated in the simplest approxima¬ 
tion where the solvent polarization is assumed to equili¬ 
brate instantaneously the propagated electronic density. 
Different schemes to account for non-equilibrium effects 
in the real-time solvent response have been published in 
the literature^^ and can be brought into our formal¬ 
ism. This will be a subject of forthcoming investigation. 
The accuracy of our real-space methodology is demon¬ 
strated by the excellent agreement of our calculations 
with the results computed by using atomic basis sets with 
the Gamess 2£ software with the advantage of allowing us 
to explore the excited states dynamics beyond the linear 
response regime with no further effort. 

The paper is organized as follows. In Sec. |TT]we explain 
the theoretical framework used through out the article. 
We start by introducing the basics of the PCM and the 
main equation to calculate the apparent surface charges 
by using the integral equation formulation of the PCM 
model. In Sec. Ill Al the regularized (singularity-free) sol¬ 
vent reaction potential in real space is derived altogether 
with the mathematical appendix [A] The PCM terms en¬ 
tering the Kohn-Sham Hamiltonian and the expression to 
calculate the free energy of the solute-solvent system are 
obtained in Sec. cm Sec. Ill Cl clarifies the approxima¬ 
tion used to calculate the solvent response in real-time. 
The numerical results calculated by using the new PCM 
implementation are discussed and benchmarked in Sec. 
IIIII Finally, in Sec. IIVI we summarize the main conclu¬ 
sions and comment on possible directions in which this 
work can be extended. 


II. THEORY 

In the framework of PCM the interaction with the sol¬ 
vent molecules is approximated by the effective Hamilto¬ 
nian, 


n = + y int , (i) 

where is the electronic Hamiltonian of the molecule 
in vacuo and K lnt is the solute-solvent interaction oper¬ 
ator expressed in terms of the reaction potentials Vfjr(r) 
and V^r) accounting for the polarization of the dielec¬ 
tric by the electronic p e { r) and nuclear p n ( r) charge den¬ 
sities of the solute system: 


N e 

Vi n t = _Y / [V^(r l ) + Vi(r l )] 

i= 1 

-i N a toms 

+ 2 E ^Vr(R/), ( 2 ) 

i=i 

where r,, refer to the electrons coordinates and R/ and 
Zj indicate the positions and the charges of the atomic 
nuclei. 

The reaction potentials in Eq. m can be defined by 
using an auxiliary apparent surface charge (ASC) a e / n (s) 
that spreads on the cavity surface T £ 



a e / n { s) 

l r — s l 


ds. 


( 3 ) 


In practice, to solve the PCM equations the continuous 
surface T needs to be discretized and approximated by 
a set of T finite surface elements commonly referred as 
tesserae. As a consequence, the integral in Eq. © trans¬ 
forms into the discrete sum: 


T/ -e/n/ \ 

V ( r ) 


T 

E 

k=l 


a e/n (s k )A k 
l r “ s fc| 


E g e/ra ( s fc) 

lr — Sfcl 
k =l 1 


( 4 ) 


where s k and A k denote, respectively, the representative 
point and the area of the fc-th tessera. In Eq. the 
area of the tesserae must be small enough to assume a 
constant value for cr e / n (s) within each element which al¬ 
lows to end up with a solvent reaction potential defined 
by the set of polarization charges {q e / n (s k )}. 

The integral equation formulation of the PCM pro¬ 
vides a general framework to calculate these polarization 
charges as: 


T 

q e/n ( s k ) = ^Q(s fc ,s;;e)y M [p e/ ”](sz), (5) 

i=i 

where Q(sfc,s/;e) denotes the PCM response matrix 
which depends on the cavity geometry and the frequency- 
dependent dielectric constant of the solvent solution. The 
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expression of the PCM matrix as well as the analyti¬ 
cal formulas to evaluate its matrix elements has been 
already published elsewhere.-’ 29 As shown in Eq. 0, 
the polarization charges by the electrons (nuclei) of the 
solute molecule are obtained by applying the PCM ma¬ 
trix to the column vector {14i[p e,/rl ](s;)} containing the 
electronic (nuclear) contribution to the molecule’s elec¬ 
trostatic potential calculated at the cavity surface as: 

V M [p e/n ]{ sz)= Jd rg^j. (6) 

In Eq. (0, p e ( r) denotes the quantum-mechanical elec¬ 
tronic density of the solute system, so that, I4i[p e ](sz) 
is nothing else than the Hartree potential at the cavity 
surface and p n ( r) = ZiS(r — Rj) is the classical 

nuclear density. 


A. The solvent reaction potential in real space 


If a real-space grid is used to represent the solvent re¬ 
action potential defined by Eq. 0 a Coulomb singular¬ 
ity will show up when the coordinate r is infinitesimally 
close to a surface representative point Sfc. This singular¬ 
ity can lead to significant numerical errors and may repre¬ 
sent an important drawback to model solvation effects in 
real-space calculations, specially for charged molecules. 
In particular, for negatively charged systems (anions), 
where the solvent response is described by a set of posi¬ 
tive polarization charges, the singular nature of the PCM 
potential in Eq. 0 might result in an artificial accumu¬ 
lation of the electrons around the cavity representative 
points which is typically manifested in a dramatic over¬ 
estimation of the solute-solvent interaction energy. The 
latter effect will be illustrated and discussed in Section 

nun 

In order to remove such a singularity, we have built a 
regularized PCM potential by approximating the set of 
point charges {<?(s*,)} by a set of Gaussian-shaped charge 
densities {p(r,Sfc)} centered at the representative points 
s*; whose widths are equal to the areas of the tesserae: 


g(r, Sfc) 


ff(Sfc) -|r- Sfc | 2 /(aA fc ) 

(tt aA k ) 3 / 2 


(7) 


In Eq. 0, g(r, s*,) is normalized to the total polarization 
charge g(sfc) and a is a numerical parameter that can be 
used to modify the width of the charge distribution. In 
the limit of a —> 0 the Eq. 0 reduces to the case of 
a point charge located at tessera representative point. 
In the appendix [0 we have calculated the electrostatic 
potential in real-space produced by the charge density 
g(r, Sfc) and obtained the following expression: 


vr(t, Sfc) 


^jU(|r 

V 7 T aAk 


Sfc| /y/aA k ), 


( 8 ) 


where Q ( x ) is a Pade approximan1>22 of order [2 /3] defined 
as: 


Q(x) 


1 +P 1 X + P 2 X 2 
1 + qix + q2X 2 + P 2 X 3 


(9) 


The coefficients in the latter equation are reported in the 
appendix El As can be noticed from Eq. @ the new 
potential u R (sfc,r) reproduces the l/|sfc — r| dependence 
for grid points located at large distances from the rep¬ 
resentative point Sfc and more important, regularizes the 
Coulomb singularity in real-space at the cavity surface 
when r = Sfc which was our primary goal (see Fig. 0. 
Finally, by summing the individual contributions of all 
tesserae we obtain an expression analogous to Eq. 0 
to calculate the regularized solvent reaction potential as 
follows: 


T 

v n(r) ='^2v R (r,Sk). ( 10 ) 

fc= l 


B. PCM terms in the real-space Kohn-Sham equations 


In the present article solvent effects in the ground- 
state electronic structure of the solute molecule are ac¬ 
counted for in the framework of Density Functional The¬ 
ory (DFT)£L. Therefore, we start from the free energy^ 
functional for the solvated system, defined as: 

G[p e ,p n ]=E°[p e } 

+ \ J d * P e ( r ) {Dt[p e ](r) + u R [p n ](r)} 

+ \f dr p n {r){v R [p e ](r)+v R [p n ](r)} , (11) 


where E°[p e ] is the total energy functional of the 
molecule in vacuo and v R (r) denotes the solvent reac¬ 
tion potential calculated in Sec. Ill Al which holds an im¬ 
plicit functional dependence on the electronic and nuclear 
charge densities throughout the polarization charges that 
are calculated following Eq. 0. If we now take the func¬ 
tional derivative on Eq. m with respect to the elec¬ 
tronic density p e { r) we obtain for the effective potential 
of the Kohn-Sham (KS) system the expression: 


v s [p e ,p n ]{r) = v 3 [p e ](r) + - {v R [p e ](r) + u R [p n ](r)} 

+ I|*'{ p -(r') + P "(r')(f^k!M, (12) 


with v ° [ p e ] (r) being the corresponding effective potential 
in vacuo including the exchange-correlation (xc) term. 
Formally, the potential u R [p e ](r) can be expressed in 
terms of the electrostatic Green functioni^, that is: 


u R [p e ](r') = |dr"G R (r',r")p e (r"). (13) 

Substituting Eq. (fl3l) in the integrand of Eq. G2D, tak¬ 
ing the functional derivative and exploiting the symmetry 
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G R (r',r") = G R (r",r') of the Green function^, the ex¬ 
pression for the effective potential v s [p e ,p n }( r) simplifies 
as follows: 

v s [p e ,p n ](v) =n s °[p e ](r) + n R [ /9 e ](r)+n R [p"](r). (14) 

The Kohn-Sham molecular orbitals (fij (r) and the eigen¬ 
values Ej are obtained by solving the following equation: 

+v s [p e ,p n ]{r) S J cpj{ r) = ejipj{ r), (15) 

The KS equations are solved numerically in real-space by 
using finite differences to approximate the kinetic energy 
operator and the derivatives of the electronic density that 
might enter the exchange-correlation potential. Compu¬ 
tational details regarding these calculations are reported 
in Sec. MB 

Eqs. (fl5l) must be solved iteratively and coupled to 
the PCM equations © and m until the equilibrium 
condition between the solute electronic density and the 
solvent polarization response has been reached. This is 
implemented by starting from a guessed electronic den¬ 
sity calculated at the level of a linear combination of 
atomic orbitals (LCAO). The latter density is plugged 
into Eq. (0 to calculate the initial approximation to the 
polarization charges {g e (sfc)}. These charges are used 
to generate the solvent reaction potential UR,(r) and af¬ 
ter solving the KS equations, the new electronic density 
is obtained from the improved orbitals. This algorithm 
is repeated so forth until self-consistency. On the other 
hand, as the nuclei coordinates remain fixed, the solvent 
response to the nuclear charge density is calculated only 
once out of the iterative scheme. 

Having solved the KS equations, the total free energy 
for the molecule in solvent can be computed as follows: 

N occ /» 

G[p e ,p n } = ~ I dr P e ( r HMP e K r ) +w R [p n ]( r )} 

3 =1 

- J dr p e (r)v xc [p e }(v) - E n [p e } + E xc [p e ] 

+ Et t [p e ,p n ), (16) 

where N occ is the number of occupied MOs, Eh is the 
Hartree energy, E xc [p e ] is exchange-correlation energy 
functional, v xc [p e ] = SE xc [p e ]/6p e ( r) and Ef^ t [p e ,p n ] is 
the contribution to the free energy due to the electro¬ 
static interaction of the solute molecule with the polar¬ 
ized solvent to be calculated from the following expres¬ 
sion: 

1 T 

E ilt\P e ,P n ) = 2 + ^M[p"](s fc )} 

k =1 

X {q e (s k ) + g”(s fc )} . (17) 


C. Solvent response in real-time 

To model the electron dynamics of the solute molecule 
under the influence of an arbitrary time-dependent ex¬ 
ternal potential, we use the time-dependent Kohn-Sham 
(TDKS) equations in the adiabatic approximation: 35 

y- +«s[p e ,p n ](r,t)^ <Pj(r,t), (18) 

where v s [p e ,p n ] is the effective potential defined by Eq. 
m evaluated in the instantaneous density p e { r, t) which 
is obtained by propagating in real-time the ground state 
KS orbitals <p(r, t = 0). 

Under these conditions, the polarization charges in¬ 
duced by the electrons of the solute molecule become also 
time-dependent and should be computed, as it is done 
for ground state calculations, coupled to the quantum- 
mechanical TDKS equations. However, in the static case 
the solvent polarization caused by the solute molecule 
is in equilibrium with the ground state density but, in 
general, as the time-dependent electronic density varies 
faster as compare with the typical timescales of the sol¬ 
vent relaxation processes, the actual polarization will lag 
behind the changing density. 

In this work, we have performed real-time TDDFT 
(RT-TDDFT) calculations for simple organic molecules 
in solvent to benchmark our real-space PCM implemen¬ 
tation. In particular, we compare the optical response 
of the solute molecules obtained from the dynamic po¬ 
larizability tensor— with linear-response TDDFT (LR- 
TDDFT) calculations performed in the frequency domain 
by using Gaussian basis sets. These calculations were 
carried out by using equal values for the dynamic and 
static dielectric constants. Under this assumption the 
calculation of the time-dependent apparent charges sim¬ 
plifies to the expression— 

T 

q e {s k ,t) = ^Q(s fe ,s ; ;e 0 )V M [p e ](si,t). (19) 

i=i 

Eq. (fl9l) is the simplest expression to calculate the sol¬ 
vent response in real-time where the dielectric polariza¬ 
tion is assumed to equilibrate instantaneously the evolved 
electronic density. This approximation, despite it ne¬ 
glects non-equilibrium effects due to the delayed compo¬ 
nents of the solvent polarization, constitutes a very good 
approach suitable to treat molecules in weakly polar sol¬ 
vents characterized in most cases by similar values of the 
static and dynamic dielectric constants. 


III. IMPLEMENTATION AND RESULTS 

We have implemented the PCM equations recalled 
above as a new module within the OCTOPUS code— In 
this section, we report on numerical calculations carried 
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out to test and benchmark the PCM implementation de¬ 
scribed in Sec. Ill At We have performed calculations of 
the ground state energies and the optical absorption re¬ 
sponse for a set of organic molecules (see Fig. [2J that 
were recently used for benchmarking purposes in Ref. 
l36l in water solution. The static and time-dependent 
KS equations m and (Hhd, have been solved in real- 
space and real-time by using OCTOPUS. Numerical de¬ 
tails about these calculations are explained in the follow¬ 
ing section. The calculated solvation free energies and 
solvatochromic shifts for selected molecules are discussed 
in sections IIIIBI and IIII Cl respectively, and compared 
with similar calculations performed with the Gamess 
software.— All the results reported in the present paper 
were computed by using the PBE generalized gradient 
approximation to the exchange and correlation energy 
functionals . 37 ’ 38 


A. Computational details 

The molecular structures shown in Fig. [2] were op¬ 
timized in gas phase by DFT calculations, as imple¬ 
mented in Gamess by using the double zeta basis set 
6-31G(d) — £2. See supplemental material at [URL will be 
inserted by AIP] for the geometries used in the reported 
calculations— 

The calculated geometries were inputed to OCTOPUS 
to perform DFT and TDDFT calculations in vacuo and 
solvent. The simulation box containing the real-space 
domain was constructed by adding spheres around each 
atomic center of radius 5 A. Calculations performed with 
larger radii showed that this value is sufficient to get to¬ 
tal and molecular orbitals energies converged for all the 
studied molecules. Ground state calculations were car- 



FIG. 2. Optimized structures of the molecules used to bench¬ 
mark the PCM implementation. The color code for the ele¬ 
ments is red for O, blue for N, grey for C, yellow for S and 
white for H. The chemical names of these molecules are re¬ 
ported in the supplemental figures S1-S4— 


ried out by using two different values for the uniform 
spacing, 0.15 A and 0.19 A, between each grid point of 
the generated mesh. Atomic pseudo-potentials have been 
used to model the core electrons for all the species. They 
were generated consistently with the Ape— code by using 
PBE approximation to the xc functional. 

To calculate the absorption spectra of some of the 
molecules shown in Fig. [2j we propagated in real-time 
the KS system from its ground state after applying, at 
t = 0, a weak delta pulse of a dipole electric field. Real 
time propagations were performed along the three polar¬ 
ization directions of the applied field by using a time step 
of At = 1.7 attosecond (as) for a total simulation time of 
20 femtoseconds (fs). The enforced time-reversal sym¬ 
metry method ) 42 ’ 43 as implemented in OCTOPUS, was 
used to propagate the time-dependent KS orbitals. The 
optical absorption strength was calculated by using the 
imaginary part of the diagonal component of the dy¬ 
namic polarizability tensor defined in the frequency do¬ 
main by Fourier transforming the time-dependent electric 
dipole. 22 The time integral over the total simulation time 
has been exponentially filtered with a damping factor of 
0.13 eV. 

The cavity surface T hosting the solute molecule is 
built as a union of interlocking spheres centered on all 
atoms but hydrogens— We have used for the radii of 
these spheres the values 2.4 A for C, 1.8 A for O, 1.9 A 
for N and 2.02 A for S— The elements of the PCM re¬ 
sponse matrix Q(sfc,sqe) are calculated by using Eqs. 
(A9, A15) and the first terms of Eqs. (A10, A14), as 
derived in the appendix of Ref. 0. We used e = 78.39 
for the static dielectric constant of water.- 4 ' Polarization 
charges are calculated with Eq. f5]) in terms of the PCM 
response matrix and the molecule’s electrostatic poten¬ 
tial at the cavity surface. The nuclear contribution to 
the latter potential can be directly evaluated by using 
Eq. ©. On the other hand, to calculate the electronic 
contribution at the cavity surface, we perform a trilinear 
interpolation by using the values of the Hartree potential 
at the grid points defining the vertices of the cube con¬ 
taining the tessera representative point (see inset in Fig. 

a. 

Once the total polarization charges are computed, the 
solvent reaction potential in real space, entering the 
Kohn-Sham equations, is calculated through Eq. m- 
All the results presented along the paper have been ob¬ 
tained by using a = 1 in Eq. (0) where the PCM poten¬ 
tial is regularized in terms of the areas of the tesserae. 
In the supplemental figures Sl-S4^ we report the weak 
dependence of the solute-solvent interaction energy for 
different values of the parameter a entering the regular¬ 
ized solvent reaction potential. 

We have used Gamess to benchmark the solvation free 
energies and the optical absorption spectra obtained with 
OCTOPUS for the investigated molecules in water. The 
same cavity surface and PCM response matrix were used 
within both codes. Total energies and excited state calcu¬ 
lations with Gamess were performed by using the triple 
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FIG. 3. Comparison of the electrostatic contribution to the 
solvation free energies calculated in real-space (RS) with Ooc- 
TOPUS with similar results obtained with Gamess by using 
Gaussian type orbitals (GTO) for the molecules shown in Fig. 
[2] The superimposed letters are used to identify each molecule 
in Table |T| 


zeta basis set 6-311+(d,p)4& 


B. Electrostatic contribution to the solvation free energy 

At the level of ground state calculations, the solvation 
free energy is an important magnitude to characterize 
quantitatively the solute-solvent interaction in thermo¬ 
dynamic equilibrium. In general, PCM allows also for 
the introduction of non-electrostatic terms^i^ in the so¬ 
lute Hamiltonian given by Eq. m which are not consid¬ 
ered in the present discussion. In what follows we focus 
our analysis on the calculated results for the electrostatic 
contribution to the solvation free energy, 

AG el = G[p s e olv , p n ] - E[pl J ac , p% (20) 

where G[p® olv ,p n ] is given by Eq. (HU) and E[p* ac , p n ] is 
the total energy of the isolated solute in vacuum. 

In Fig. [U real-space calculations for AG e i are com¬ 
pared with the corresponding results obtained by using a 
Gaussian basis set for the whole set of molecules shown 
in Fig. [2] in water solution. These results were obtained 
by using a grid spacing of 0.19 A and a = 1 to gener¬ 
ate the regularized solvent reaction potential (Eq. (flOl) ) in 
real space. From the figure, an excellent agreement be¬ 
tween the different calculations is evident. As a matter of 
fact, all plotted circles hardly deviate from the diagonal, 
showing maximum differences of about 0.3 kcal/mol. The 
calculated free energies are also reported in Table H] for 
the sake of completeness. We have performed the same 
calculations without regularizing the PCM potential by 
using Eq. fU to evaluate the solvent reaction field in real- 
space and obtained very similar results as it is shown in 
Fig. S5 of the supplemental material^ This suggests 


TABLE I. Electrostatic contribution to the solvation free en¬ 
ergies in kcal/mol plotted in Fig. [3] Real-space (RS) calcu¬ 
lations were performed by using a grid spacing of 0.19A with 
Octopus. The analogous results by using Gaussian type or¬ 
bitals (GTO) were calculated with Gamess. 



molecule 

ag2 s 

AGg TO 

a) 

mol-140 

-0.79 

-0.83 

b) 

mol-016 

-1.76 

-1.88 

c) 

mol-163 

-2.60 

-2.68 

d) 

mol-142 

-3.04 

-3.29 

e) 

mol-079 

-3.14 

-3.46 

f) 

mol-036 

-4.17 

-4.19 

g) 

mol-160 

-4.00 

-4.20 

h) 

mol-206 

-5.05 

-5.20 

i) 

mol-200 

-5.56 

-5.51 

j) 

mol-090 

-6.07 

-6.19 

k) 

mol-013 

-8.05 

-8.24 

1) 

mol-117 

-10.48 

-10.49 

m) 

mol-069 

-10.58 

-10.77 


that numerical problems due to the Coulomb singularity 
are unlikely and also that the value a = 1 is a very good 
choice for the regularization. The same analysis has been 
carried out by using a smaller grid spacing of 0.15 A and 
we found, again, the same trend and accuracy. 

However, there may be a threshold for the minimum 
distance between a tessera and a grid point for which nu¬ 
merical instabilities should appear. In order to illustrate 
this possibility, we consider one of the worst case sce¬ 
nario. In Fig. [I] we plot the interaction energy given by 
Eq. (fT71) calculated for the chlorine anion in water as a 
function of the distance d m ; n between the representative 
point s at the cavity surface and a grid point iq. This 



FIG. 4. Calculated electrostatic interaction energy for a chlo¬ 
rine anion in water for different values of the distance between 
the tessera representative point s and the real-space point ri . 
0.15 A was used as grid spacing. The F surface in this case is 
defined by a sphere of radius 2.17A. 
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is sketched by the inset superimposed to Fig. [4] where 
{]"!, ..,r 8 } label the vertices of the cube surrounding the 
tessera. Results are presented for two cases: open circles 
correspond to calculations performed with a PCM po¬ 
tential given by Eq. 0; solid circles plot the interaction 
energies computed by using the regularized potential as 
defined in Eq. (USD- As it appears from the plot, there is 
a critical distance (in this case d m i n ~ 2.8 x 10 -4 A) be¬ 
low which the singular nature of the Coulomb potential 
manifests dramatically. However, as it can be no¬ 
ticed from the same figure, such numerical problems are 
completely absent in the calculations carried out with the 
regularized PCM potential predicting accurate energies 
for all distances, even when the tessera coordinates fully 
overlaps with the grid point (d m i n = 0). It is noteworthy 
that this result was also obtained by setting a = 1. 


C. Solvent effects in the optical absorption response 

In this section we report and discuss the optical 
response of selected molecules: piperazine (mol-013), 
trimethylaminne (mol-016), biphenyl (mol-079) and ni¬ 
trobenzene (mol-206). The absorption spectra of these 



3 3.5 4 4.5 5 5.5 6 6.5 

energy (eV) 

FIG. 5. Absorption spectra of nitrobenzene molecule (mol- 
206) calculated with (a) RT-TDDFT by using OCTOPUS and 
(b) LR-TDDFT with Gamess. 


systems have been calculated in gas phase and in water 
by solving in real-time and real-space the TDKS (Eq. 
(ITHll 'l and PCM (Eq. (fliil) ) equations simultaneously, as 
implemented in our code of choice, OCTOPUS. The com¬ 
puted spectra have been benchmarked with LR-TDDFT 
calculations by using Gaussian type orbitals as imple¬ 
mented in Gamess. 

In Fig. [5] we plot the calculated strength function ver¬ 
sus the excitation energy for the nitrobenzene molecule. 
The spectra obtained by doing RT-TDDFT calculations 
are shown in the upper panel while results obtained 
within linear response approximation are reported in the 
lower panel. Notice from this figure that the calculated 
spectra, both in vacuo and in solvent, with RT-TDDFT 
and LR-TDDFT are almost identical, that is, the main 
absorption bands appears at the same excitation ener¬ 
gies with similar relative intensities. In particular, we 
focus our analysis on how the solvatochromic shifts of 
the main absorption bands obtained by real-time prop¬ 
agation compare with the analogous result obtained by 
using a Gaussian basis set; this comparison offers a di¬ 
rect measure of the quality of the PCM implementation 
presented along this paper. To that aim, we have la¬ 
beled with roman numerals the absorption peaks in the 
gas phase spectrum that can be clearly identified in the 
solvated one for the calculated molecules. See supple¬ 
mental material at [URL will be inserted by AIP] where 
the absorption spectra of the molecules piperazine (mol- 
013), trimethylaminne (mol-016) and byphcnil (mol-079) 
are plotted in Figs. S6-S8»^ 

In Fig. E we correlate the solvatochromic shifts for 
the main absorption bands of the investigated molecules 
as obtained from RT-TDDFT calculations with the cor¬ 
responding values calculated with Gamess in the linear 
response approach. The plotted values report for the en- 



FIG. 6. Solvatochromic (solv.) shifts of the main absorption 
bands of the molecules mol-013, mol-016, mol-079 and mol- 
206 in water. Real-time (RT) TDDFT calculations are corre¬ 
lated with analogous results obtained in the linear response 
(LR) approximation. See also Table [III 





















TABLE II. Solvatochromic shifts in eV of the main absorption 
bands for the selected molecules in water obtained from real¬ 
time (RT) and linear response (LR) TDDFT calculations. 


molecule 

peak 

solv. shift RT 

solv. shift LR 

mol-013 

I 

0.22 

0.23 


II 

0.02 

0.03 


III 

0.09 

0.09 


IV 

0.04 

0.07 

mol-016 

I 

0.10 

0.11 


II 

-0.10 

-0.10 


III 

0.17 

0.15 


IV 

0.08 

0.07 

mol-079 

I 

-0.19 

-0.19 


II 

-0.03 

-0.03 


III 

-0.16 

-0.16 

mol-206 

I 

-0.43 

-0.44 


II 

-0.17 

-0.18 


III 

-0.28 

-0.28 


ergy differences taken at the absorption maxima between 
the solvated and gas phase spectra of the calculated 
molecules. As it turns out from this plot, RT-TDDFT 
calculations predict shifted excitation energies due to the 
dynamical solvent polarization in perfect agreement with 
linear response calculations performed by using Gaussian 
atomic functions. We have found a maximum deviation 
of 0.03 eV from the diagonal which corresponds to the 
solvatochromic shift of the highest-energy band of piper¬ 
azine. For the sake of completeness we report in Table ITIl 
the values of the solvatochromic shifts plotted in Fig. [G] 

IV. CONCLUSIONS 

In this work we presented a new methodology to ac¬ 
count for solvation effects in the electronic and optical 
properties of molecular systems by using density func¬ 
tional approximations in real-space and real-time. We 
have used the Integral Equation Formalism of the Po¬ 
larizable Continuum Model to calculate the apparent 
charges at the cavity surface hosting the solute molecule. 
To prevent numerical instabilities in real-space calcula¬ 
tions due to Coulomb singularities in the solvent reaction 
potential we introduced a set of Gaussian functions, cen¬ 
tered at the tesserae representative points, to distribute 
the polarization charges within the area of each surface 
element. This allowed us to derive a simple analyti¬ 
cal expression to evaluate the PCM potential over the 
whole simulation domain. We used the regularized sol¬ 
vent potential to do ground state calculations by solv¬ 
ing the self-consistent Kohn-Sham (KS) equations by us¬ 
ing the PBE generalized gradient approximation to the 
xc energy. The new PCM implementation was applied 
to compute in real-space the solvation free energies of a 
set of organic molecules in water by using the OCTOPUS 


code. For benchmarking proposes we compared with sim¬ 
ilar calculations by using Gaussian atomic orbitals per¬ 
formed with Gamess and found an excellent agreement 
between both results that showed maximum differences 
of about 0.3 kcal/mol. 

We have also extended the PCM equations to the real¬ 
time domain to investigate solvent effects in the electron 
dynamics under the influence of a time-dependent exter¬ 
nal field. We calculated the optical spectra of selected 
molecules in water by propagating in real-time the KS 
orbitals and the PCM apparent charges, after perturb¬ 
ing the molecule with a weak instantaneous dipole. The 
obtained solvatochromic shifts for the main absorption 
bands compare extremely well with the analogous results 
calculated in the frequency domain within the linear- 
response approximation by using Gamess. At present, 
the real-time solvation effects has been included by as¬ 
suming that solvent polarization equilibrates instanta¬ 
neously the evolved electronic density. This approach, 
despite it neglects non-equilibrium effects in the dielec¬ 
tric response, is a very good approximation to investigate 
molecules in weakly polar solvents. 

The methodological developments and the results pre¬ 
sented in this paper set the grounds for further exten¬ 
sions to model solvation effects in the framework of real- 
space and real-time electronic structure calculations by 
using density functional approximations. For instance, 
the inclusion of non-equilibrium effects to describe the 
solvent response in real-time can be incorporated into 
our computational machinery as proposed in Ref. [23]. 
This would contribute to evaluate dielectric effects more 
rigorously in the time-resolved optical spectroscopy of 
molecular systems in solution to probe ultrafast charge 
transfer and non-linear optical phenomena. Moreover, 
this work constitutes a good starting point to add solva¬ 
tion effects to model the coupled electron-nuclei dynam¬ 
ics in the Ehrenfest approximation as implemented in the 
Octopus code4&i2 
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Appendix A: Regularization of the solvent reaction 
potential in real-space 

In real-space calculations, it is evident from Eq. © 
that the solvent reaction potential Wi,(r) shows a singu¬ 
larity for grid points that are infinitesimally close to the 
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solute cavity surface, that is, for r —> s k . To regularize 
the PCM potential at this limit, we distributed isotropi¬ 
cally the polarization charges within each tessera by using 
a Gaussian function centered at the representative point 
Si and width equals to the area A k , 

g( r, s fc ) = Ne -W-s k \ 2 /(<*A k )' ( A1 ) 


The integrand of the latter equation can be further sim¬ 
plified by changing to the dimensionless variable x 2 = 
y 2 /(aAk) and scaling out the factor ait to the integra¬ 
tion limit, that is: 


v(r, Sfc) = T(\r - s k \/y/aA k ), (A9) 

\/naAk 


The parameter a in Eq. CU) has been introduced to 
modify the width of the Gaussian smearing out the po¬ 
larization charges. The constant N can be obtained from 
the normalization condition: 


J dr p(r, Sfc) = A^47 t J dx x 2 e ^ kX = q( Sfc), 


(A2) 


with the function /F(x) defined as, 


T(x) 



T(3/2) — T(3/2,a:' 2 ) | 
x' 2 J ' 


(A10) 


The argument of function J-{x) takes values between 


where /3k = 1 /(aAk) and x = |r — s* |. Making use of the 
parametric integral, 


dx e 


— @kX 2 


1 j TV 

2 V /3fc ’ 


we obtain the normalized charge density: 


p(r,Sfc) 


g( s fc) p — |r—Sfc | 2 / (aA k ) 

(tt aA k ) 3 / 2 


(A3) 


(A4) 


On the other hand, by using the Gauss law and the 
spherical symmetry of g , the radial component of the 
electric field at a certain distance x k = |r — Sfc| from the 
representative point, can be calculated as follows: 



E (xk) = j (A5) 
47T60 Sfc 

where q(x k ) is the amount of polarization charge inside 
the spherical surface of radius x k to be calculated by 
integrating the charge density in Eq. tiHD , 

l * dv ^’ Vi ‘ AkK < A6 > 

By changing to the variable t = y 2 /(aAk) the latter in¬ 
tegral rewrittes as follows: 

q(x k ) = r dt t^ 2 e-* - f°° dt ^/ 2 e-* 

v 7r •'0 Jx\l{aA k ) 

= ~^ [r(3/2) - r(3/2, x 2 /{aA k ))\ , (A7) 

v 7r 


FIG. 7. Calculated values for the function J-(x) (solid circles) 
obtained from Eq. (1 A10I) . The Pade approximant G{x), used 
to fit the asymptotic behavior and functional dependence of 
3F(x ), is plotted as a solid line. 

0, for the case r = Sfc and oo. In Fig. [T] we plot the 
calculated values for the integral in Eq. (IA10I) as func¬ 
tion of x in the interval between 0 and 10 atomic units. 
Notice that the potential resulting from Eq. (IA9I) re¬ 
produce the expected 1 jr dependence at large distances 
from the tesserae representative points, while in the limit 
of r = Sfc regularize the Coulomb singularity to the con¬ 
stant value 2q(&k)/y/naAk which was our primary goal. 
These asymptotic limits as well as the functional depen¬ 
dence of the F(x) can be perfectly fitted by using the 
Pade approximant >211 


where T(s) and T(s,a;) denote the complete and upper 
incomplete gamma functions^, respectively. Now, with 
an analytic expression for q(x k ) we can calculate the elec¬ 
trostatic potential generated by g(r,s k ) at any point in 
real space by integrating Eq. (EH), 


u (r, Sfc) 


2g(sfc) 



dy 


f T(3/2) — T(3/2, y 2 /(aAk)) 

l y 2 


■ (AS) 


G{x) 


1 + PlX + P2X 2 

1 + qix + qix 2 + P 2 X 3 ’ 


(All) 


with the coefficients pi = 0.119763, P 2 = 0.205117, q\ = 
0.137546 and 52 = 0.434344. The excellent agreement 
between the proposed fit and the exact results calculated 
from the integral in Eq. (IA10I) is shown in Fig. 0 
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